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ABSTRACT 

We investigate the clustering of dark matter halos in Lagrangian space in terms of their two-point 
correlation function, spanning more than 4 orders of magnitudes of halo masses. Analyzing a set of 
collisionless scale-free 128 3 -particle N-body simulations with spectral indices n = —2,-1, we measure 
the first two Lagrangian bias parameters b\ and 62 relating halo and mass correlations. We find that 
the Mo & White leading-order formula for b\ describes the clustering of halos with mass M ;> M* quite 
accurately, where M* indicates the characteristic non-linear mass. Smaller halos turn out to be less 
clustered in Lagrangian space than predicted by Mo & White. Our findings are consistent with the 
recent results of Jing for the clustering of halo populations in Eulerian space, demonstrating that the 
discrepancies between the A-body and analytical Mo & White prediction for the bias exist already in 
Lagrangian space. This shows that a more refined theoretical algorithm for selecting halos in the initial 
conditions needs to be developed. Finally, we present a very accurate fitting formula for the linear halo 
bias factor b\ in Lagrangian space. 

Subject headings: cosmology: theory - galaxies: statistics - large-scale structure of the Universe 



1. INTRODUCTION 

Virialized dark matter halos in the Universe are not 
distributed in the same way as the underlying dark mat- 
ter. This is true whether one takes the final positions 
of the halos in Eulerian space or their initial positions in 
Lagrangian space (Mo & White 1996 (MW); Catelan et 
al. 1998 (CLMP); Catelan, Matarrese & Porciani 1998; 
Jing 1998; Sheth & Lemson 1998). The numerical tour 
de force by Jing (1998), who thoroughly investigated with 
unmatched accuracy the clustering of dark matter halos in 
Eulerian space, demonstrated explicitly that i) the halo-to- 
mass bias is independent of the halo separation (at least in 
the scale-free case n — — 2 and in the linear regime); ii) the 
MW Eulerian linear bias correctly describes the clustering 
of halos of masses M^M+, but systematically underpre- 
dicts it for any value of the spectral index n if M <; M+, 
where M* is the typical non-linear mass. However, it is 
impossible to understand solely on the basis of the Eule- 
rian investigation whether the discrepancies between the 
numerical results and the analytical MW predictions are 
due to a failure of the algorithm for identifying the halo po- 
sitions in Lagrangian space, or to the effects of non-linear 
shear dynamics (not accounted for in the original MW ap- 
proach) on the mapping of halo positions from Lagrangian 
to Eulerian space, or to a combination of the two. 

In this Letter, we employ scale-free collisionless A-body 
simulations (i.e. with density parameter Q — 1 and initial 
power spectra P{k) cx k n ) to investigate the clustering of 
dark matter halos in Lagrangian space in terms of their 
two-point correlation function, spanning more than 4 or- 
ders of magnitudes in halo masses. We compare the halo 
correlation function to the correlation of the underlying 
dark matter for n = — 2, — 1. 

We find that the theoretical 'underclustering' reported 
by Jing for masses M <J M* is already present in the initial 



conditions as well (but as 'overclustering', since the first- 
order bias is negative for small masses), and cannot be 
due exclusively to the subsequent non-linear effects of the 
shear dynamics acting on small scales and unaccounted for 
in the original spherical collapse model of MW. We then 
argue that the standard Press-Schechter approach ('ex- 
tended' or not), on which Mo & White (1996) and Cate- 
lan et al. (1998) based their speculations, is inadequate 
for identifying the locations of halos in the initial condi- 
tions, and a major effort should be devoted to fin- ding an 
improved algorithm. § 2 reports the details of the present 
investigation. Finally, an accurate fitting formula for the 
bias b\ in Lagrangian space is given, which should be con- 
sidered as the Lagrangian version of Jing's fitting formula 
for the Eulerian case. §3 contains our conclusions. 

2. HALO CLUSTERING FROM JV-BODY SIMULATIONS 

2.1. Simulations and Halos 

The simulations used here are similar to those of Lacey 
& Cole (1994), who used them to test halo merging histo- 
ries. They were performed using the AP 3 M code of Couch- 
man (1991) with 128 3 particles. The force-softening used 
was L/1280, with L the size of the periodic box. Initial 
positions and velocities were generated by displacing par- 
ticles from a uniform 128 3 grid according to the Zel'dovich 
approximation, assuming an initial scale-free power spec- 
trum and Gaussian statistics. We considered 4 realizations 
for 2 different spectral indices n = —2,-1. For each simu- 
lation we recorded, for many epochs, positions and veloci- 
ties of all particles. The output times were chosen so that 
increased by a factor \/2 between subsequent output 
times. 

For each output time we selected dark matter halos in 
the simulations employing the 'friends-of-friends' group 
finder with a linking length equal to 20% of the mean 
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interparticle distance (e.g. Davis et al. 1985). We checked 
that results obtained using the spherical overdensity group 
finder (Lacey & Cole 1994) are essentially identical. We 
excluded halos containing less than 20 particles or more 
than 20000 particles. We moved all the particles belonging 
to a given halo back to their initial (Lagrangian) positions, 
then computed the position of their centre of mass. We 
used the latter as the 'halo position' in Lagrangian space. 
In such a way, for each output time, we constructed a 
catalogue of halos indicating their mass and position in 
Lagrangian space. 

2.2. Halo correlation function in Lagrangian space 

We computed the mean correlation function ^ between 
halos in a given mass interval, where a bar denotes mass 
averaged quantities. Self-similar scaling allowed us to com- 
bine data from different output times in order to reduce 
the Poisson fluctuations due to the finite number of halos 
within the box. We considered every output time contai- 
ning more than 100 halos in the same mass interval, and 
we binned the distributions of halo separations in units of 
r/R+, where R% oc M*. Finally, we computed the halo 
correlation function using the estimator 

i i 

where the index i runs over different output times Zi of 
the same simulation, Ni(r/R+) is the number of halo pairs 
in the i-ih. output, and Nf° 1 (r/R ir ) is the corresponding 
quantity for a Poisson process with the same number den- 
sity. The average symbol ( ) is introduced since we con- 
sidered information coming from the different temporal 
outputs. In this way, for each mass interval, we collected 
four realizations of the Lagrangian correlation function of 
dark matter halos. 

This procedure allows us to achieve two goals: i) to re- 
duce statistical fluctuations by increasing the number of 
halo pairs; ii) to extend the mass interval and the range 
of halo separations that may be sampled. In fact, the box 
size L and the minimum halo mass in the simulation are 
fixed, while i?* and M* increase with time. This means 
that the correlation function for halos with M ^ M+ is 
measured mainly from the early output times, while that 
for halos less massive than M* comes mostly from later 
output times. In order to simulate bootstrap resampling, 
we assigned as the standard error the Poisson errorbar 
multiplied by a factor V3 (Mo, Jing & Borner 1992). 

2.3. Lagrangian bias parameters 

Now we want to test whether the Lagrangian halo cor- 
relation is related that of the mass (£ m ) through a re- 
lation (see Catelan et al. 1998; Porciani et al. 1998 and 
references therein) 




where the symbol bi denotes the i-th Lagrangian bias 
factor, and £ m is calculated according to linear theory. 
For instance, the 'extended' Press-Schechter approach 
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Fig. 1. — Dots: Best fit values of Sj (top panel) and b?, (bottom 
panel) as a function of variance and halo mass on the same scale 
for n = —2. Errorbars denote one-sigma uncertainties. Long-dashed 
line: MW and CLMP Lagrangian bias with <5 C = 1.686. Short- 
dashed line: Jing's fit to the Eulerian bias minus 1. Continuous 
line: our fitting formula for b\ given in the main text. Both the 
theoretical predictions and the fitting formula are averaged over the 
Press-Schechter mass function. Note how the halo correlation re- 
mains non-zero at M = M* through the effect of 62 , even though 61 
vanishes there. The horizontal errorbar in the bottom-right corner 
of top panel indicates the mass intervals considered. The box inside 
the top panel shows the best fit correction f(o\. ) to the MW formula 
for 61, when M < M*: log/ = (0.34 ± 0.11) - (0.73 ± 0.12) log cr^ . 

(Bond et al. 1991) leads to the MW and CLMP ex- 
pressions for b\ and &2, namely b\ = S c /a^ — 1/S C and 

62 = (Sc/^m ~ ^)/°m' wnere °m indicates the mass vari- 
ance on scale M (Cole & Kaiser 1989; Mo & White 1996; 
Mo, Jing & White 1997). Note that, by definition, a M = S c 
for M = M*. To test this biasing model against simula- 
tions, we have to consider a finite range of halo masses. 
Eq.(j|) then implies £ h = b\£ m + p^m + . . ., where b k is 
the mean of b^ in the mass interval, weighted by the mass 
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2. — As in Fig. 1 but for n = —1. In this case, the best fit for 
log / = (0.13 ± 0.13) - (0.41 ± 0.10) logCT^ 



function. In order to obtain statistically reliable values for 
b\ and 62, we compute the averages (£ m ) and directly 
from the initial conditions of the simulation, following the 
same averaging procedure we used for the halo correlations 
(in this case it corresponds to averaging over output times 
using the volume of the r/R+ bin as a weight). Note that, 
after averaging, generally differs from (£ m ) 2 , so that 
eq.(§) implies (&) = 6f(£ ro ) + |b|<^) + . . .. This avera- 
ging procedure is important to account for the finiteness 
of the box: in this way, the lack of any Fourier components 
of the density would be equally experienced by both mass 
and halo distributions. Finally, to check the reliability of 
eq.(Q) we apply the x 2 test, assuming that every is 
normally distributed around a mean value. 

Halos are objects of finite size and we expect exclusion 
effects to dominate £/, at separations of the order of the 
halo Lagrangian size Rl- In a sample composed of identi- 
cal spherical halos, spatial exclusiveness implies £^ = — 1 
for r < 2Rl and a strong compensating positive correla- 
tion at r 2Rl- The magnitude of this positive corre- 



lation is related to the volume fraction occupied by the 
spheres, since J °° drr 2 £fc(r) = 0. Inspired by this, we de- 
cided to exclude from the x 2 test all data corresponding to 
separations r < 2i? max , where R max = i?*(M max /il4) 1 / 3 is 
the characteristic Lagrangian radius of the largest halo in 
the mass interval considered; the exact value of the mini- 
mum separation considered has been determined on a case 
by case basis, checking for the stability of our results with 
increasing the number of data points at small separations. 
We avoided averaging between different realizations of the 
same power spectrum: the entire data set of up to 4 values 
of for each r/R+ interval is considered in the x 2 test- 
Minimization of \ 2 with non-holonomic constraints b\ > 
and b\ > is used to obtain estimates of the bias pa- 
rameters. The resulting Xmin shows that the a priori bias 
hypothesis-the bias relation in eq.(^)-is acceptable to a 
95% confidence level over a wide interval of halo masses. 
However, for n = —2, samples with M > 16M* are in- 
consistent with eq. (j2j) at the 99% confidence level. This 
is caused by the presence of noticeable differences (espe- 
cially at very large separations) between the halo corre- 
lation functions extracted from different realizations. In- 
deed, very massive halos are exponentially rare and more 
affected by statistical fluctuations. On the other hand, we 
could also have underestimated the uncertainty in com- 
puting (&). 

In our approach a small x 2 is found only if eq.pft with 
the same values for b\ and 62 accurately describes every 
clustering realization, for a given n. We accounted for 
both the dispersion between realizations (since data from 
all the realizations enter the x 2 test ) and the uncertainty 
of the single data points (since the residuals between fit 
and data are divided by the simulated bootstrap errors). 
This approach is different and much more severe than first 
averaging £/, and £ m over the realizations and finally fit- 
ting the data using a least-squares method. In this case, a 
larger dispersion between realizations creates larger error- 
bars for the averaged quantities, thus making x 2 smaller. 
In Jing (1998) no statistical test is performed: the bias 
relation at large separations, where £ m <C 1, is assumed to 
be linear. For each separation and realization, the squared 
Eulerian bias is estimated by computing £h/£m- These va- 
lues are then averaged over different realizations and the 
estimated dispersion of the mean value is taken as the er- 
rorbar. Finally, by considering different separations, it is 
checked whether the averaged biases have consistent va- 
lues, within their statistical uncertainties. 

Figures 1 and 2 show the best-fit results for b\ and b\. 
The errorbars are obtained by projecting along the 6-axes 
the contours of constant x 2 at Xmin + 1- Concerning b\, the 
comparison against theoretical predictions shows a quanti- 
tatively good agreement for M > M+. However, the square 
bias parameter of lower mass halos is significantly smaller 
than predicted by Mo & White (meaning that the bias b\ 
is less negative than MW value). Findings of this kind 
have been published by Jing (1998), who investigated the 
Eulerian halo clustering in the quasi-linear regime, giving 
an accurate fitting formula for the Eulerian analogue of 61 . 
Our result shows incontrovertably that the effects discove- 
red by Jing are already present in Lagrangian space. 

There is no rigorous way to apply Jing's bias fitting 
formula in Lagrangian space, since the non-local dyna- 
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mics of the mass density field enters the transformation. 
However, as a first approach, and only to check the order 
of magnitude of the effect, we can follow Mo & White 
(1996) in assuming a spherical evolution of the coarse- 
grained mass density held, to map from Lagrangian to 
Eulcrian space. In this case, the first Eulerian bias term 
is given by 1 + b\. Jing's formula, b\ = bj mg — 1, turns 
out to be a reasonably good description of our data, sho- 
wing that the 'underclustering' of small halos in Eulerian 
space detected by Jing is correspondingly already present 
in their Lagrangian clustering; the intervening dynamics 
is approximately well accounted for by the standard map- 
ping from Lagrangian to Eulerian space in laminar regime. 
We elaborate below our own fitting formula for the linear 
Lagrangian bias, which turns out to be accurate to 10% 
over the entire mass range investigated (with the exception 
of the zero crossing region for b\). 

Our results for 62 are of course less conclusive than 
for bi, mainly because of the larger uncertainties. How- 
ever, they confirm the information extracted considering 
bi . Quantitatively, the CLMP expression for b^ gives good 
estimate for halos with M > M+, while it overpredicts the 
numerical outcome for b\ when M < M*. For M ~ M*, 
our results are a factor of ~ 1.4 smaller than predicted 
by the CLMP formula. This is the range in which our 
method is best suited to compute 62 , since b\ vanishes and 
the halo correlation function is proportional to No 
firm conclusion can be drawn for n = — 1, since practically 
the_entire data set only allows one to set an upper limit 
on b\. It is not unlikely that, for the less negative spectral 
indices, bigger simulations are better suited to quantify 
the halo clustering up to second-order biasing, above all 
in Lagrangian space. It will be worthwhile to address this 
point in a future work. 

2.4. Accurate fitting formula for b\ 

Our results for b\ can be accurately parametrized by in- 
troducing a mass-dependent, multiplicative correction to 
the MW formula, namely 

bl(cr M )=f(a M )(5 c /**-l/5 c ) 2 . (3) 

For n — —2, the original Mo & White formula (i.e. / = 1) 
with 5 C — 1.686, as suggested by the spherical collapse 
model, describes the data for M > M* to within 10% ac- 
curacy. Smaller masses instead require / = 2.19/c^- 46 . 
For n = — 1 and M > M* our numerical data are very well 
described by the MW formula for b\ with a lower collapse 
threshold S c ~ 1.52. In this case, for M < we obtain 
/ = 1.35/cr^ 82 . These fitting formulae are extremely ac- 
curate in describing our data set. Their simple power-law 
behavior encourages further theoretical investigation. For 
instance, with data for more than two spectral indices n, 
one may attempt to fit the dependence of b\ on n. 



3. DISCUSSION AND CONCLUSIONS 

Employing 128 3 -body scale-free simulations, we ana- 
lyzed the clustering of dark matter halos in Lagrangian 
space. The main results of this investigation can be sum- 
marized as follows: i) assuming a correlation model as 
in eq.(2), the first two Lagrangian bias factors b\ and 62 
are strongly mass-dependent over the 4 orders of magni- 
tude in mass investigated; it) the clustering of halos with 
mass above the non-linear mass, M ;> M*, is fairly well de- 
scribed by the MW formula for the linear Lagrangian bias, 
both for n = — 1 and n = —2; in) halos with non-linear 
masses M <; M* are less clustered (have a smaller corre- 
lation amplitude) than what the leading order Lagrangian 
bias of Mo & White would predict. 

When these results are combined with the ones recently 
obtained by Jing (1998) about the clustering of halos in 
Eulerian space, we can disentangle the question of whether 
the discrepancies between Jing's numerical results and the 
Mo & White theoretical predictions are mainly due to 
the effects of the non-linear shear dynamics, effective on 
smaller scales, or to a possible failure of the halo selection 
algorithm in the initial conditions-a question actually left 
unsolved by Jing. Clearly, since, as we showed, the same 
effects discovered by Jing are essentially already present 
in Lagrangian space, our investigation suggests that it is 
time to improve on the classical Press-Schechter algorithm 
for identifying halos in Lagrangian space: since it assumes 
spherically symmetric collapse, it is not surprising that it 
fails in correctly counting the small halo masses, where 
departures from spherical collapse can be cosmologically 
relevant. The failure of the Mo & White formula for the 
Lagrangian and Eulerian bias is presumably related to 
the departure of the halo mass function from the Press- 
Schechter form at low masses, e.g. Lacey & Cole (1994). 

Finally, we derived a fitting formula for the linear La- 
grangian bias b\ that is relevant for accurately predicting 
the clustering of dark matter halos, above all in the low- 
mass tail. Our fitting formula can be considered as the 
Lagrangian equivalent of eq.(3) in Jing (1998). These re- 
sults are highly relevant for predicting the clustering of low 
luminosity galaxies, most of which lie in lower-mass halos 
(e.g. Baugh ct al. 1998). Modelling of galaxy clustering 
found in present and future galaxy redshift surveys will 
provide an important application of the present results. 
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